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Chapter  1 


Introduction:  source  of  turbulent  flow 
inhomogeneities  in  the  oceanic  flows 


1  •  1  Introduction 

The  passage  of  a  coherent  electromagnetic  beam  through  a  turbulent  medium  results 
in  a  change  of  light  velocity  which  in  turn  causes  distortion  in  intensity  and  phase 
of  the  beam.  Light  propagation  through  the  Earth’s  atmosphere  has  been  studied 
extensively  and  many  models  have  been  developed  to  describe  atmospheric  turbu¬ 
lence  (see  for  example  the  review  by  Strohbehn  [67]).  To  date  similar  studies  of  the 
propagation  of  light  in  turbulent  water  are  very  sparse  [71].  Field  measurements 
of  small-angle  scattering  [55]  show  that  the  volume  scattering  function  exhibits  a 
sharp  peak  which  is  orders  of  magnitude  greater  than  that  obtained  from  either  the 
laboratory  measurements  of  Spinrad  [63]  or  Mie  type  calculations  for  non-turbulent 
conditions  [62].  This  effect  (Le.  peaking  of  volume  scattering  function  at  small 
angles)  has  been  attributed  to  scattering  by  turbulence-induced  inhomogeneities  in 
the  refractive  index  of  sea  water  [77].  One  of  the  main  goals  of  this  work  is  to  verify 
this  conjecture. 

Currently,  increased  use  of  optical  detectors  in  the  marine  environment  (for  exam¬ 
ple  biological  oceanography  and  in  underwater  imaging  and  communications)  makes 
quantifying  the  role  of  turbulence  an  important  tcisk.  Fortunately,  only  small-scale 
fluctuations  in  the  refractive  index  due  to  the  flow  field  are  relevant  for  light  propa¬ 
gation.  These  small-scale  features  are  only  weakly  dependent  on  the  large-scale  flow 
and  thus  are  similar  for  different  flow  types  [3]. 
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The  interaction  of  light  with  fluid  inhomogeneities  associated  with  turbulence 
has  long  been  utilized  in  laboratory  observations  of  turbulence  (t.  e.,  using  the  shad¬ 
owgraph  technique)  [73].  This  technique  has  also  been  used  in  situ,  e.g.  to  visualize 
a  double-diffusive  instability  [42]. 

It  is  extremely  difficult  to  accurately  measure  nearforward  light  scattering  due 
to  the  existence  of  the  primary  light  beam;  likewise  any  simplifying  assumptions 
about  the  structure  of  turbulent  flow  fields  can  lead  to  unexpected  errors  [72].  Thus 
I  have  chosen  to  use  numerical  simulation  to  investigate  light  propagation  through 
a  turbulent  medium. 


1.2  Inhomogeneities  in  the  refractive  index 

Inhomogeneities  in  the  refractive  index  (IRI)  are  characteristic  of  turbulent  flows. 
Under  typical  oceanic  conditions,  these  inhomogeneities  appear  to  be  temperature 
dominated.  In  general,  however,  the  real  part  of  the  refractive  index  of  seawater 
varies  with  changes  of  temperature,  salinity,  and  pressure  (changes  in  the  imaginary 
part  of  the  refractive  index  of  water  are  negligibly  small  compared  with  the  real 
part  [8]).  The  effect  of  pressure  can  be  neglected  for  spatial  scales  on  the  order  of  a 
meter  or  smaller.  The  variance  of  the  IRI  (assuming  that  salinity  and  temperature 
fields  are  not  correlated)  can  be  expressed  as  the  sum  of  the  variance  associated  with 
temperature  and  salinity:  (n^)  ^  =  {n^)  -|-  (n|)^^^,  where  ()  denotes  spatial  aver- 

^  is  the  refractive  index,  and  nj  and  ns  are  the  contributions  of  temperature 
and  salinity  respectively. 

In  the  ocean,  background  temperature  and  salinity  induce  quite  different  variance 
in  the  IRI.  For  example,  analysis  of  data  from  the  coastal  region  off  Oregon  (Figure 
3.5)  show  that  variance  of  the  refractive  index  due  to  temperature  is  much  larger 
than  that  due  to  salinity  (J.  Mourn,  pers.  comm.).  Thus,  only  temperature-induced 
IRI  will  be  considered  in  this  study,  so  hereafter  n  =  ny.  The  contribution  of  the  IRI 
to  scattering  has  been  postulated  to  be  largest  for  the  smallest  IRI  structures  [77]. 
Therefore,  in  this  study  I  will  concentrate  on  light  scattering  by  IRI  structures  at 
the  smallest  scale  of  temperature  field.  Since  the  behavior  of  temperature  embedded 
in  the  turbulent  field  can  be  numerically  simulated  and  then  exhaustively  studied. 
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we  start  with  direct  numerical  simulation  (DNS)  of  the  temperature  field  in  the 
turbulent  flow. 
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Chapter  2 


Direct  numerical  simulations  of  passive  scalars 
with  Pr  >  1  advected  by  turbulent  flow 


2.1  Introduction 

The  spectral  behavior  of  kinetic  energy  spectra  in  homogeneous,  isotropic  turbu¬ 
lence  is  fairly  well  established  through  results  from  experiments,  theoretical  studies, 
and  numerical  simulations.  These  studies  have  been  reviewed  in  several  monographs 
[49,  43].  For  flows  at  high  Reynolds  numbers  the  most  important  paxt  of  the  spec¬ 
trum  is  the  k  inertial  subrange,  which  according  to  the  Kolmogorov  1941  theory 
terminates  below  the  Kolmogorov  or  dissipation  wavenumber 

kK  =  {elu^yi*  =  (2.1) 

where  e  is  the  kinetic  energy  dissipation  rate,  u  is  the  kinematic  viscosity,  and  tjk  is 
the  Kolmogorov  length  scale.  Experimental  results  indicate  that  the  inertial  range 
ends  around  k  (0.1  —  0.2)/:^.  Many  theoretical  and  numerical  results  [39,  37, 
18,  12,  44]  support  the  following  form  for  the  dissipation  range  for  wavenumbers 
k  >  0.2kj( 

E{k) {kfkK)°exp{-ak/kK).  (2.2) 

The  exponential  is  well  established  and  has  been  confirmed  by  the  experiments  of 
Sreenivasan  [65]  but  the  algebraic  prefactor  is  uncertain  though  several  investigations 
suggest  a  «  —2  in  the  range  0.2  <  k/k^  <  4.  While  the  inertial  subrange  exists  only 
for  sufficiently  high  Reynolds  number  turbulence,  generally  inaccessible  to  numerical 
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simulations,  the  dissipation  range  form  (2.2)  appears  to  be  universally  valid  for  all 
Reynolds  numbers  and  can  be  easily  simulated  numerically.  The  universality  in 
particular  implies  that  those  features  of  turbulence  that  are  determined  by  scales 
from  the  vicinity  of  kK  can  be  investigated  through  direct  numerical  simulations 
(DNS)  irrespective  of  Reynolds  number.  This  observation  is  used  later  in  our  work. 

In  many  physical  situations  one  encounters  passive  scalaxs  advected  by  turbulent 
velocity  fields  like  temperature  and  pollutants  in  air  and  water.  Following  Kol¬ 
mogorov  s  local  isotropy  hypothesis  Obukhov  [54]  and  Corrsin  [14]  predicted  that 
scalar  variance  spectra  at  high  Reynolds  numbers  will  also  have  a  inertial  sub¬ 
range.  The  behavior  of  the  variance  spectra  for  wavenumbers  beyond  the  Corrsin- 
Obukhov  range  depends  on  the  Prandtl  (or  Schmidt)  number  Pr  =  iz/zc,  where  k 
is  the  scalar  diffusivity.  For  Pr  <\  Batchelor  et  al.  [4]  predict  the  following  form 
for  the  spectrum  of  the  variance  of  a  scalar  advected  by  the  velocity  scales  from  the 
inertial  subrange 

Ee{k)  ~  (2.3) 

where  x  is  the  scalar  variance  dissipation  rate;  the  formula  is  valid  for  wavenumbers 
k  greater  than  the  Obukhov-Corrsin  wavenumber,  koc 

koc  =  =  Pr^/'kK.  (2.4) 

This  form  of  the  spectrum  has  been  confirmed  by  the  numerical  simulations  of  Chas- 
nov  et.  al.  [11], 

In  oceanographic  applications  temperature  and  salinity  are  important  scalars, 
and  both  have  Pr  >  1.  The  case  of  Pr  1  was  investigated  theoretically  by 
Batchelor  [3]  who  derived  the  following  form  of  the  scalar  spectrum  in  the  viscous- 
convective  range 

Eg{k)  = -^k~^  exp{Kk'^f'y),  (2.5) 

where  7  is  the  average  value  of  the  least  principal  rate  of  strain  (a  negative  quantity) 
acting  on  scales  k  >  ki{ .  The  derivation  is  based  on  a  physical  picture  of  the  scalar 
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field  at  scales  k  ^  kfc  being  deformed  by  an  essentially  uniform  gradient  of  the 
velocity  at  much  larger  scales  k  ^  kx-  For  k  less  than  the  Batchelor  wavenumber 

ka  =  {cluK^yi*  =  Pr^l'^kK  =  1/j7b,  (2.6) 

the  exponential  factor  in  (2.5)  is  approximately  equal  to  unity  and  the  algebraic 
prefactor  k  ^  dominates  such  that  the  so-called  Batchelor  spectrum  becomes 

Ee{k)  =  ~k-\  (2.7) 

The  final  form  of  the  Batchelor  spectrum  is  determined  by  specifying  the  strain  rate 
7.  For  flows  at  high  Reynolds  numbers  it  is  assumed  that  a  good  estimate  is  given 
by  the  rate-of-strain  of  Kolmogorov  eddies 

7  =  --(c/i/)^/^,  (2.8) 

q 

where  ^  is  a  universal  constant.  Using  available  experimental  data  Batchelor  [3] 
estimated  g  «  2.  Gibson  [22]  provides  a  range  of  \^ues  \/3  <  q  <  2-\/3  deduced 
from  the  incompressibility  and  kinematics  of  fluid  elements.  Oceanic  measurements 
provided  values  of  9  =  3.9  ±  1.5  [27]  and  q  =  3.7±  1.5  [53]. 

Each  of  the  previous  spectral  predictions  for  scalar  spectra  results  from  a  sepa¬ 
rate  set  of  physical  assumptions  making  the  predictions  independent  of  each  other. 
Gibson  [21]  and  [22]  proposed  a  unified  theory  of  passive  scalars  advected  by  turbu¬ 
lence  valid  for  all  values  of  the  Prandtl  number.  The  theory  is  based  on  the  physical 
mechanism  of  generation  of  maxima  and  minima  of  the  scalar  concentration  by  the 
local  rates  of  strain.  This  mechanism  is  distinctly  different  from  the  mechanism 
invoked  by  Batchelor  for  Pr  I .  While  for  Pr  >  1  Gibson’s  prediction  for  the 
spectra  of  the  scalar  is  the  same  as  Batchelor’s  prediction  (2.5),  for  Pr  1  the 
Gibson’s  theory  implies  the  intermediate  range  for  koc  <k<  ka  before 
range  develops  for  k  >  ka-  The  evidence  in  support  of  Gibson’s  theory  is  discussed 
by  Clay  [13],  and  Kerr  [37].  In  the  context  of  our  work  Gibson’s  analysis  has  two 
appealing  features.  First,  it  suggests  that  the  Batchelor  k~^  range  may  appear  even 
if  the  Prandtl  number  is  not  much  greater  than  unity  (i.e.  the  condition  Pr  >  1 
replaces  Batchelor’s  Pr  ^  1).  Second,  it  implies  a  universal  scaling  independent  of 
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the  Prandtl  number  for  scalar  spectra  in  the  scalar  dissipation  range.  This  so-called 
Batchelor  scaling  is  obtained  by  normalizing  wavenumbers  by  the  Batchelor  length 

scale  rjB  and  spectra  by  Under  this  scaling  expression  (2.5)  is  given  by 

a  self-similar  form 


q{kr]B)  *  exp(-g(fc7;B)2). 


(2.9) 


Kraichnan  [40]  considered  the  effects  of  fluctuations  of  the  rates  of  strain  in  space 
and  time  in  Batchelor’s  analysis  and  showed  that  while  k~^  behavior  survives,  the 
Gaussian  factor  in  (2.5)  is  not  universal  and  is  replaced  by  a  simple  exponential.  He 
also  obtained  the  same  result  using  Lagrangian-History-Direct-Interaction  (LDHI) 
closure  theory.  Kraichnan ’s  equation  for  a  scalar  spectrum  was  solved  in  a  closed 
form  by  Mjolsness  [48]  giving  the  following  result 


E0{k)  =  5^fc  ^(1  -f  ^)  exp(-^),  (2.10) 

where  ^  =  (30/c/A)^/^A;,  and  the  constant  A  can  be  in  principle  calculated  from  the 
theory.  By  comparing  functional  expressions  (2.10)  and  (2.7)  in  the  limit  A:  -4  0  I 
get  A  =  (5/q)(c/i/)^/^.  Using  the  Batchelor  scaling  formula  (2.10)  becomes 


EejkrjB) 


—  q{kr}B)  ‘(1  +  i^qY^^kriB) exp{-{6qY^^kT]B). 


(2.11) 


Similar  expressions,  containing  k~^  range  and  the  exponential  form,  are  also  pre¬ 
dicted  by  other  analytical  theories  of  turbulence  [52,  29,  58].  The  analytical  theories 
also  make  predictions  about  the  universal  constant  q  in  (2.11).  Kraichnan  [40]  es¬ 
timates  q  <  0.9  using  LHDIA,  Newman  and  Herring  [52]  give  q  =  1.68  using  Test 
Field  Model  and  Qian  [58]  gives  q  =  2-\/E.  Only  the  last  value  is  consistent  with  the 
experimental  evidence. 

In  summary  the  classical  picture  of  the  scalar  spectrum  in  high  Reynolds  number 
flow  and  Pr  ^  1  as  uniformly  presented  in  turbulence  literature  [70,  49,  43]  is  the 
inertial  k  behavior  for  k  <  the  Batchelor  k~^  behavior  for  kn  <  k  <  kB,  and 
the  exponential  decay  for  k  >  kB- 
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Despite  the  wealth  of  theoretical  results  in  favor  of  this  picture  the  experimen¬ 
tal  evidence  is  ambiguous,  ranging  from  supporting  to  rejecting  in  part  or  entirely 
the  spectral  forms,  in  particular  the  k~^  Batchelor  range,  arrived  at  by  theoretical 
consideration.  Gibson  and  Schwarz  [24]  observed  the  Batchelor  spectrum  for  tem¬ 
perature  and  salinity  in  laboratory  measurements  in  water  and  the  approximate 
behavior  for  temperature  spectrum  is  also  suggested  by  oceanographic  field  measure¬ 
ments  of  Grant  et  al.  [27]  and  Oakey  [53].  However,  later  measurements  of  Gargett 
[20]  contradicted  the  universality  of  the  Batchelor  spectrum  and  cast  doubt  on  ear¬ 
lier  experimental  results.  Gibson  [23]  attributes  Gargett’s  negative  results  to  her 
failure  to  account  for  large  intermittence  of  oceanic  turbulence.  The  oceanographic 
implications  of  the  discrepancy  are  summarized  in  the  review  of  Gregg  [28]. 

More  recently  Miller  and  Dimotakis  [47]  reported  results  of  laboratory  measure¬ 
ments  of  scalar  spectra  for  high  Schmidt  number  (~  1.9  •  10^)  in  a  turbulent  jet 
which  does  not  show  any  evidence  of  the  k~^  spectral  behavior.  Note  that  according 
to  the  Batchelor  theory  the  large  value  of  Schmidt  number  1)  in  this  experiment 
should  allow  for  an  easy  detection  of  the  k~^  range.  The  k~^  range  is  also  absent  in 
the  high  Schmidt  number  2-D  mixing  experiments  of  Gollub  et  al.  [25]. 

In  addition,  when  the  Batchelor  range  is  in  fact  observed  the  quality  of  the 
experimental  data  often  does  not  allow  a  clear  distinction  between  the  spectral  forms 
(2.11)  and  (2.9)  in  the  fax-dissipation  range.  Consequently  experimental  results 
usually  are  interpreted  assuming  the  earlier  derived  Batchelor  spectrum  (2.9)  and 
no  comparison  studies  with  the  Kraichnan  spectrum  are  available.  The  practical 
importance  of  these  spectral  expressions  lies  in  the  fact  that  all  scalar  fluctuations 
and  scalar  dissipation  are  effectively  determined  by  scales  from  the  k~^  Batchelor 
range.  The  dissipation  rates  in  turn  determine  the  mixing  coefficients  for  scalars  [17] 
which  are  critical  to  understand  the  small-scale  physics  of  the  oceans  as  well  as  the 
large  scale  circulation  and  global  climate.  The  knowledge  of  spatial  power  spectra  of 
temperature  fluctuations  at  small  scales  is  also  needed  to  address  problems  of  sound 
and  light  propagation  in  water  [69].  Several  scalar  models  have  been  examined  for 
those  problems  [30]  but  again  lack  of  experimental  data  at  the  largest  wavenumbers 
introduces  uncertainties  in  the  analyses.  In  view  of  the  practical  importance  of 
the  Batchelor  spectrum  and  the  scatter  in  the  predicted  values  of  the  universal 
constant  q  as  well  as  experimental  controversies  concerning  the  existence  of  the  the 
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k  '  range  itself,  further  work  in  this  area  is  clearly  needed.  Carefully  designed 
numerical  simulations  may  complement  experiments  and  theory  in  contributing  to 
the  resolution  of  these  controversies. 

While  numerical  simulations  have  been  used  extensively  to  investigate  scalar 
spectra  for  Pr  <  I  [37,  61,  57]  no  similar  numerical  effort  has  been  devoted  to  the 
case  of  Pr  >  1.  The  numerically  discouraging  aspect  of  this  case  is  that  substan¬ 
tially  higher  resolution  is  needed  to  simulate  the  scalar  field  than  the  velocity  field. 
Indeed,  the  theoretically  suggested  location  of  the  Batchelor  spectrum,  k  >  kK,  im¬ 
plies  that  the  simulations  should  extend  far  into  the  dissipation  range  to  values  of 
k':^  kK-  This  effectively  limits  the  turbulence  simulations  on  current  supercomput¬ 
ers  to  very  low  values  of  the  Reynolds  number.  However,  two  observations  suggest 
that  it  may  be  possible  to  obtain  the  Batchelor  spectrum  in  a  properly  designed  nu¬ 
merical  experiment  despite  the  inability  of  DNS  to  simulate  high  Reynolds  number 
flows.  First,  as  already  argued  by  Batchelor  [3]  the  form  of  the  relation  (2.5)  should 
also  hold  for  low  Reynolds  number  flows,  as  long  as  Pr  >  1  and  there  exist  large 
velocity  scales  providing  uniform  rates  of  strain.  More  generally,  the  existence  of 
the  Batchelor  range  for  scalars  with  Pr  >  1  can  be  viewed  as  a  consequence  of  the 
scalar  diffusion  in  a  strain  field  that  is  spatially  fairly  uniform  and  random.  This  is 
persuasively  demonstrated  by  the  simple  physical  model  of  Antonsen  et  al.  [2]  where 
the  extensive  k  ^  range  was  obtained  by  numerically  solving  the  passive  scalar  equa¬ 
tion  with  the  velocity  field  prescribed  as  a  long  monochromatic  wave  with  a  random 
phase.  A  similar  approach  has  been  employed  by  Holzer  and  Siggia  [33]  for  two- 
dimensional  flows  and  by  Pumir  [57]  for  three-dimensional  flows.  In  both  cases  only 
a  few  low  wavenumber  velocity  modes,  i.e.  large  scales,  actively  advect  the  passive 
scalar  which  is  allowed  to  develop  scales  much  smaller  than  the  forcing  velocity  field. 
Holzer  and  Siggia  [33]  observed  the  k~^  range  for  simulations  performed  with  512^ 
modes  and  hyperdiffusion  while  Pumir  [57]  estimates  that  in  three  dimensions  the 
resolution  of  1000^  modes  would  be  needed  to  observe  the  k~^  range.  Also  Metals 
and  Lesieur  [46]  attribute  an  ‘anomalous’  passive  scalar  range  in  their  large  eddy 
simulations  to  stirring  of  temperature  by  the  large  energy-containing  eddies.  While 
in  many  of  these  cases  the  k  ^  range  is  observed,  the  prevailing  rate  of  strain  7  in 
(2.7)  is  generally  different  from  the  estimate  (2.8),  expected  to  hold  only  for  high 
Reynolds  number  turbulence.  Second,  as  noted  above,  there  is  evidence  that  the 
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energy  dissipation  range  (2.2)  is  independent  of  (or  only  weakly  dependent  on)  the 
Reynolds  number.  Therefore,  the  spectral  dynamics  of  scales  from  the  vicinity  of 
the  Kolmogorov  wavenumber  kx ,  where  the  dissipation  range  is  located,  should  be 
the  same  for  large  and  small  Reynolds  number  flows.  When  Kolmogorov  scaling 
is  used  the  increase  in  the  Reynolds  number  only  extends  the  energy  spectrum  to 
smaller  wavenumbers  without  affecting  scales  from  the  vicinity  of  the  Kolmogorov 
scale.  Since  Batchelor  theory  assumes  that  the  scalar  variance  dynamics  is  driven 
by  the  Kolmogorov  velocity  scales,  the  universality  of  the  energy  dissipation  spectra 
implies  that  the  velocity  field  for  a  low  Reynolds  number  turbulent  flow  should  have 
the  same  effect  on  the  scalar  in  the  Batchelor  range  as  a  high  Reynolds  number 
flow.  This  observation  suggests  that  it  is  plausible  to  expect  not  only  the  same 
k  ‘  spectral  behavior,  but  also  the  same  scaling,  in  particular  for  7,  to  hold  for  all 
Reynolds  numbers.  Note,  however,  that  previous  numerical  simulations  performed 
for  this  problem  [33,  57,  2]  purposely  did  not  attempt  to  simulate  the  dissipation 
range  of  the  velocity  field  and  thus  were  unable  to  address  this  issue. 

In  view  of  these  arguments  I  believe  that  it  is  worthwhile  to  explore  the  dynamics 
of  passive  scalars  for  Pr  >  1  in  low  Reynolds  number  turbulent  flows  by  numerically 
solving  the  full  Navier-Stokes  equation  for  velocity  and  a  passive  scalar. 

2.2  Basic  Equations  and  Numerical  Methods 

The  flow  dynamics  are  described  by  the  incompressible  Navier-Stokes  equations  for 
velocity  u(x,  t)  and  the  transport  equation  for  the  passive  scalar  ^(x,  i) 

dn  1 

■^  +  (u- V)u= --Vp-f-i/V^u-|-F„,  (2.12) 

^  •  u  =  0,  (2.13) 

Qj.  +  (u  •  V)0  =  kV  9  -H  F$.  (2.14) 

In  the  equations  ^  is  a  constant  density,  p  is  the  pressure,  and  1/  and  k  are  the 
kinematic  viscosity  and  the  molecular  diffusivity,  respectively.  The  terms  and 

Fe  signify  forcing  added  to  the  equations  in  order  to  obtain  statistically  stationary 
turbulence. 
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e  flow  ,s  assumed  to  be  contained  in  a  cube  of  side  L  =  2n  and  periodic 

boundary  condrt.ons  in  all  three  spatial  directions  are  imposed  on  the  independent 

vanables.  The  domain  is  discretised  in  physical  space  using  N  uniformly  spaced  grid 

points  m  each  direction  resulting  in  a  mesh  size  Ai  =  L/N  and  a  total  of  grid 

points.  The  equations  of  motion  are  transformed  to  spectral  space  using  the  discrete 
Fourier  transform 

E  exp(-ik  •  x),  (2.15) 

and  the  inverse  transform  is 


u(xj  =  ^u(k)exp(fk-x), 

k 


(2.16) 


where  x  are  the  mesh  points  in  physical  space  and  k  are  the  discrete  wavenumbers 
with  components  h,  =  ±„,Ai,  n,  =  0, 1. 2,  ...iV/2,  i  =  1,2,3,  and  =  2x/i  =  I, 
The  distinction  between  physical  and  spectral  representation  for  a  given  quantity  is 
made  through  its  argument  x  or  k,  respectively. 

The  equations  are  solved  using  a  pseudo-spectral  numerical  method  of  Rogallo 
[60]  in  the  implementation  of  Yeung  and  Pope  [76].  I  have  employed  the  forcing 
scheme  of  Sullivan  et  al.  [68]  in  which  the  total  energy  of  several  low  wavenumber 
modes  IS  kept  constant  while  evolution  of  individual  modes  through  nonlinear  inter¬ 
actions,  subject  to  the  global  constraint,  is  allowed.  Specifically,  the  sum  of  squared 
amplitudes  of  modes  in  a  sphere  of  radius  Kj  =  2.5AA:  is  kept  constant  for  both  the 
velocity  and  the  passive  scalar.  This  is  accomplished  by  multiplying  all  modes  in 
the  forced  sphere  by  the  same  constant  factor,  usually  not  larger  than  1.02,  at  the 

end  of  each  time  step.  This  restores  the  energy  (or  the  scalar  variance)  in  the  sphere 
to  the  value  at  the  beginning  of  the  time  step. 

Physical  quantities  of  interest  for  isotropic  turbulence  are  described  in  terms  of 

the  scalar  wavenumber  k  =  |k|  by  averaging  over  thin  spherical  shells  defined  for  an 
arbitrary  quantity  /(k)  as 

<W>  =  lE/(k).  (2.17) 

K  \f 
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where  (...)  denotes  the  shell  average  and  the  summation  extends  over  all  modes 
in  the  shell  of  thickness  Ak  centered  at  A:  =  |k|.  The  energy  and  the  passive  scalar 
spectra  are  defined  as  follows 

E{k)  =  4rf(i!i„(k)<(k)),  (2.18) 


£«(i)  =  4ri^{<l(k)r(k)),  (2.19) 

and  the  corresponding  dissipation  spectra  are 

D{k)  =  2uk^E{k\  (2.20) 


Dg{k)  =  2Kk'^Ee{k).  (2.21) 

The  integral  of  E{k)  over  k  gives  turbulent  kinetic  energy  per  unit  mass  3/2^2  = 
3/2«  ,  where  the  overbar  denotes  averaging  in  physical  space  and  u'  is  the  rms 
turbulent  velocity,  and  the  integral  of  Ee{k)  gives  the  scalar  variance  ¥.  The  inte¬ 
grated  dissipation  spectra  give  the  dissipation  rates  of  the  kinetic  energy,  e,  and  of 
the  scalar,  x-  The  Taylor  microscale  is  computed  as  A  =  (15u'V/e)^/2 
croscale  Reynolds  number  is  Rx  =  u'X/ u.  An  important  time  scale  for  the  evolution 

of  turbulence  is  the  large  eddy  turnover  time  T,  =  where  Tp  is  the  integral 

length  scale 


(2-22) 

Lp  is  used  to  define  the  macroscale  Reynolds  number  Re  =  u'Lpju  and  the  Peclet 
number  Pe  =  u'Lp/K  =  Pr  Re. 


Values  of  the  velocity-derivative  skewness  and  the  mixed-derivative  skewness 
Sue  are  often  used  in  direct  numerical  simulations  to  assess  how  well  turbulence  is 
developed.  The  velocity  derivative  skewness  Su  is  given  by  the  formula 


(2.23) 
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For  isotropic  turbulence  the  vorticity  production  rate  is  proportional  to  the  velocity- 
derivative  skewness  The  mixed-derivative  skewness 


(fe)  (&)^) 


(2.24) 


IS  related  to  the  nonlinear  transfer  of  scalar  variance  to  small  scales  and  takes  a  value 
zero  when  there  is  no  net  cascade  to  higher  wavenumbers.  For  statistically  steady 

isotropic  turbulence  and  may  be  expressed  entirely  in  terms  of  energy  and 
passive  scalar  spectra  [36] 


{e/l5i/)2 

(2.25) 

fk^m)dk 

15  (e/15r/)j(x/6«)' 

(2.26) 

2.3  Results 

Calculations  were  carried  out  on  the  Cray  C90  with  resolution  of  up  to  240^  modes. 
In  addition  to  the  turbulent  velocity  field  in  each  run  we  carried  three  passive 
scalars  with  values  of  the  Prandtl  number  of  3,  5,  and  7.  The  low  wavenumber 
modes  for  both  velocity  and  scalars  were  forced  as  described  earlier.  As  an  initial 
condition  I  used  the  prescribed  energy  spectrum  proportional  to  the  exponential 
E[k,0')  C'exp(  ck),  where  constants  C  and  c  were  chosen  for  a  given  resolution 
to  provide  a  prescribed  value  of  the  microscale  Reynolds  number.  This  form  of  the 
initial  conditions  is  consistent  with  the  experimental  decay  of  the  dissipation  spectra 
at  larger  wavenumbers  (2.2).  Since  low  wavenumber  modes  Q  <  k  <  Kj  —  2.5 Ak 
are  forced  and  excluded  from  the  analysis,  I  did  not  attempt  to  match  in  the  ini¬ 
tial  conditions  the  algebraic  prefactor  in  (2.2).  Equivalently,  I  could  have  initialized 
simulations  with  a  constant  energy  spectrum  for  the  forced  modes  k  <  K/  and 
vanishing  spectrum  for  k  >  Kj,  but  such  simulations  would  require  longer  time  to 
reach  a  statistically  stationary  state.  The  initial  scalar  spectra  were  proportional 
to  the  energy  spectrum.  The  simulations  were  started  with  the  resolution  of  128^ 
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modes  and  run  for  about  one  large  eddy  turnover  time.  Subsequently,  the  gen¬ 
erated  dataset  was  used  as  a  restart  for  a  run  at  higher  resolution  162^  which  was 
allowed  to  evolve  until  a  steady  state  for  the  scalar  field  was  achieved  (it  takes  longer 
for  the  scalar  spectra  to  achieve  equilibrium  than  for  the  energy  spectra  [57]).  In 
this  procedure  at  the  restart  time,  the  high  wavenumber  modes,  outside  the  initial 
(smaller)  spectral  domain,  are  all  zero.  After  the  restart  these  modes  evolve  rapidly 
as  they  gam  energy  from  the  modes  in  the  original  smaller  spectral  domain  by  the 
nonlinear  interactions.  The  change  in  the  Reynolds  number  was  accomplished  by 
changing  the  viscosity.  I  found  that  this  approach  in  increasing  resolution  minimizes 
the  CPU  time  needed  to  generate  a  fully  developed,  statistically  steady  spectrum 
at  the  target  resolution  while  still  generating  proper  results.  Turbulence  parameters 
in  all  runs  were  chosen  such  that  all  energy  and  scalar  spectra  axe  well  resolved  in 
the  final  state.  Due  to  computer  time  constraints  the  highest  resolution  case  was 
run  for  about  one  large  eddy  turnover  time  after  restarting  from  one  of  the  162^ 
runs.  Towards  the  end  of  this  run  the  high  wavenumber  modes,  initially  zero  at 
the  restart  time,  still  experience  slight  energy  increase  indicating  imbalance  between 
the  nonlinear  transfer  and  the  dissipation  in  this  part  of  the  spectrum.  However, 
the  remaining  parts  of  the  spectrum  are  found  in  the  statistical  equilibrium  and  the 
velocity  derivative  skewness  and  the  mixed-derivative  skewness  for  this  run  attain 
generally  accepted  equilibrium  values.  These  observations  imply  that  simulated  tur¬ 
bulence  is  sufficiently  well  developed  for  the  purpose  of  our  analysis  which  focuses  on 
the  behavior  at  low  wavenumbers.  Only  results  from  the  simulations  performed 
with  the  higher  resolutions  of  162^  and  240^  modes  are  reported  here.  Parameters 
for  these  runs  are  gathered  in  Table  1  and  Table  2. 
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Run 

Resolution 

1/ 

k 

^max 

wm 

A 

VK 

■■ 

e 

Rex 

Re 

xl62 

0.02 

76 

1.55 

~  1.1 

mi 

0.052 

~  25 

34 

0.60 

aal62 

0.01 

76 

1.34 

~  0.8 

|I9 

0.047 

~36 

62 

0.46 

a240 

0.0033 

113 

1.09 

~  0.4 

■wWl 

0.080 

~  77 

185 

0.51 

Table  2.1:  Parameters  of  the  simulations  for  the  velocity  field. 


Run 

Pr 

VB 

X 

Pe 

■«IK» 

102 

mg 

xl62 

■iliBil 

170 

in 

7 

0.042 

QQ 

238 

0.45 

3 

igg 

186 

aal62 

5 

■wikM 

310 

mi 

7 

0.026 

434 

0.38 

3 

0.26 

555 

0.41 

a240 

5 

0.26 

925 

0.39 

7 

0.0098 

0.27 

1295 

0.36 

Table  2.2:  Parameters  of  the  simulations  for  the  scalar  field. 

Results  for  the  spectra  presented  below  are  calculated  using  shell  averages  (2.17) 
with  the  first  two  forced  shells  removed  and  are  suitably  nondimensionalized  bringing 
data  from  all  runs  to  a  self-similar  form.  The  classical  Kolmogorov  scaling  proce¬ 
dure  for  the  energy  spectra  [37]  uses  'qK  and  for  nondimensionalization  of 

wavenumbers  and  spectra,  respectively.  For  scalar  spectra  the  Batchelor  [3]  scal¬ 
ing  is  usually  employed  where  qs  and  are  used  to  nondimensionalize 

wavenumber  k  and  the  spectra,  respectively. 

2.3.1  The  velocity  spectra 

Fig.  2.1  presents  log-log  plot  of  three-dimensional  energy  spectra  using  the  Kol¬ 
mogorov  scaling  and  multiplied  by  {kqi^yf^.  For  wavenumbers  kqfc  <  0.2  the  plotted 
functions  appear  flat  which  suggests  a  5/3  slope  for  the  energy  spectra  characteristic 
of  the  inertial  range.  This  range  of  small  wavenumbers  is  enlarged  in  Fig.  2.2  where 
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the  vertical  axis  is  plotted  using  a  linear  scale.  The  approximate  inertial  range  is 
clearly  present  in  run  a240  where  the  Kolmogorov  constant  Ck  is  between  1.2  and 
1.6  for  krjK  <  0.2.  This  value  is  consistent  with  the  values  observed  experimen¬ 
tally  but  smaller  than  reported  in  numerical  simulations  where  values  around  2.0 
are  frequently  found.  No  inertial  range  behavior  is  found  for  runs  at  lower  Reynolds 
numbers.  It  is  also  clear  that  the  Kolmogorov  scaling  succeeds  in  bringing  spectra 
at  different  Reynolds  numbers  into  a  self-similar  form.  The  values  of  the  velocity 
derivative  skewness  for  all  runs  are  given  in  Table  1  and  are  consistent  with 
the  generally  accepted  value  of  -0.5  [36]  for  well  developed  turbulence.  Fig.  2.3 
presents  the  same  quantity  as  in  the  previous  figures  but  plotted  on  the  log-linear 
scale  which  emphasizes  the  far  dissipation  range.  In  the  fax  dissipation  range  the 
energy  spectrum  is  believed  to  be  a  simple  exponential,  i.e  E{k)  ~  B  exp{—aT]Kk). 
Our  data  give  value  of  B  =  8.5  and  a  =  5.4  consistent  with  values  of  B  =  8.4  ±  0.6 
and  G  =  5.1  ±  0.1  of  Kida  and  Murakami  [38].  This  fit  is  particularly  good  for  the 
coefficient  B  due  to  large  resolution  of  our  runs  and  relatively  poorer  for  the  slope 
a  because  of  the  early  termination  of  the  a240  run. 

In  Fig.  2.4  I  plot  the  normalized  kinetic  energy  dissipation  spectra.  In  all  cases 
the  dissipation  is  concentrated  at  krjK  around  (1  -  2)  x  10"^  For  the  highest  reso¬ 
lution  case  I  observe  the  well  defined  dissipation  peak  but  for  the  two  cases  at  the 
lower  resolution  the  maximum  of  the  dissipation  spectrum  occurs  at  the  smallest 
unforced  wavenumbers.  A  complete  overlap  between  the  energy  containing  and  the 
dissipation  range  is  consistent  with  the  lack  of  the  inertial  subrange  in  the  latter 
cases. 

2.3.2  The  scalar  variance  spectra 

In  the  previous  subsection  I  showed  that  the  energy  spectra  are  self-similar  and  are 
consistent  with  other  numerical  and  experimental  results.  This  suggests  that  the 
numerically  simulated  velocity  fields  can  be  used  to  reliably  advect  passive  scalars 
in  order  to  study  their  dynamics.  This  suggestion  is  supported  by  the  observed 
values  of  the  mixed-derivative  skewness  in  the  simulations.  One  of  the  secondary 
predictions  of  Gibson’s  theory  [22]  is  that  the  mixed  derivative  skewness  is  a  constant 
independent  of  the  Reynolds  and  Prandtl  numbers.  Numerical  simulations  of  Kerr 
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Figure  2.1:  Normalized  energy  spectra  E{krjK){kTiKfl^lUu^Yl^  for 
aal62,  o-run  a240,  *-run  xl62). 


all  runs  (x-run 
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Figure  2.2:  Normalized  energy  spectra  E{kT]K){kr]Kfl^ !{  in  the  range  of  low 

wavenumbers.  The  horizontal  lines  in  the  plot  mark  the  range  of  values  of  the 
Kolmogorov  constant  (x-run  aal62,  o-run  a240,  *-run  xl62). 
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Figure  2.3:  Normalized  energy  spectra  E[krif^)[krif^)^^^ f plotted  using  log- 
linear  scale  to  emphasize  the  exponential  form  (x-run  aal62,  o-run  a240,  *-run 
xl62). 
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Figure  2.4:  The  normalized  energy  dissipation  spectra  E{kr]K){kT}KY in  the 
range  of  low  wavenumbers  (x-run  aal62,  o-run  a240,  *-run  xl62)  . 
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[36,  37]  agree  with  this  prediction  giving  S^e  —  —0.5.  The  mixed  derivative  skewness 
is  found  to  be  around  -0.4  in  the  experiments  of  Clay  [13].  In  our  simulations  Sus  is 
close  to  -0.4  for  all  analyzed  runs  (see  Table  2).  This  also  implies  that  the  interactions 
between  the  velocity  and  scalar  are  well  developed  in  the  final  states  analyzed  in  this 
paper. 

The  controversial  evidence  for  the  existence  of  the  Batchelor  spectrum  comes 
mostly  from  few  laboratory  and  oceanographic  measurements  and  numerical  solu¬ 
tions  of  simplified  model  problems  [10,  53,  24,  33,  57,  2].  To  our  knowledge  there  are 
no  existing  numerical  simulations  of  full  Navier-Stokes  equations  for  passive  scalar 
fields  with  Pr  >  1  which  exhibit  the  Batchelor  range  behavior. 

In  Fig.  2.5  I  plot  the  scalar  spectra  multiplied  by  krjB  using  the  Batchelor 
scaling  for  all  three  runs  and  all  three  Prandtl  numbers  Pr  =  3, 5, 7.  All  curves 
collapse  tightly  to  one  universal  form  which  shows  that  the  Batchelor  scaling  is 
indeed  universal  for  Pr  >  1  in  the  dissipation  and  the  far  dissipation  range. 

Because  of  the  multiplication  by  kr]s  the  flat  part  of  the  curves  in  the  range  of 
small  wavenumbers  0.02  <  krig  <  0.2  corresponds  to  the  Batchelor  range  k~^  be¬ 
havior  which  is  observed  in  our  simulations  even  for  scalar  with  the  Prandtl  number 
as  low  as  3. 

In  Fig.  2.6  I  plot  the  normalized  dissipation  spectra  Es(kriB){krjB)^ I 
for  all  cases.  The  Batchelor  scaling  clearly  collapses  the  dissipation  spectra  for  all 
Reynolds  numbers  and  Prandtl  numbers  and  the  location  of  the  dissipation  peak 
around  kTjB  =  0.25  is  consistent  with  oceanic  data  of  Oakey  [53]. 


2.3.3  Parameters  of  the  scalar  spectra  in  the  dissipation 
range 

In  the  spectral  forms  proposed  by  Batchelor  (2.9)  and  by  Kraichnan  (2.11)  the  con¬ 
stant  q  is  unknown  while  all  other  parameters  are  either  prescribed  (i/  and  k)  or  can 
be  computed  from  the  data  (e  and  x).  Values  of  the  universal  constant  q  are  obtained 
from  the  least  square  fit  of  formulas  (2.9)  and  (2.11)  to  the  numerical  data  for  the 
run  a240  giving  qe  =  3.9  ±  0.25  and  =  5.26  ±  0.25,  respectively.  In  Fig.  2.7  I 
present  the  results  of  the  fit  compared  with  the  numerical  data.  Kraichnan’s  spectral 
form  describes  our  data  extremely  well,  with  some  exception  in  the  exponential  tail 
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corresponding  to  the  scalar  with  Pr  =  7.  This  can  be  attributed  to  the  previously 
discussed  fact  that  in  this  case  the  scalar  field  still  evolves  in  the  range  of  the  largest 
wavenumbers.  In  the  case  of  the  Batchelor  form  the  best  fit  curve  Fig.  (2.7)  shows 
large  systematic  deviation  from  our  numerical  data.  The  systematic  deviations  of 
the  Batchelor  prediction  from  experimental  data  measured  in  situ  has  long  been 
noted  in  oceanographic  literature.  Oakey  [53]  observed  that  in  the  k~^  portion  of 
the  spectrum  the  best  fit  to  the  Batchelor  spectrum  systematically  underpredicts 
the  experimental  data.  I  observe  similar  behavior  in  our  simulations  (see  Fig.  2.7), 
with  the  numerical  data  located  above  the  Batchelor  curve  for  wavenumbers  before 
the  dissipation  peak.  Our  estimate  of  qb  =  3.9  is  in  a  good  agreement  with  the  value 
9  ~  ^  obtained  by  Oakey  [53]  from  oceanic  measurements. 

Sometimes  in  oceanographic  experiments  the  kinetic  energy  dissipation  rate  c 
is  determined  indirectly  from  the  measured  temperature  variance  spectra.  This  is 
accomplished  by  assuming  that  the  temperature  spectra  are  properly  described  by 
the  Batchelor  formula  (2.9)  with  the  assumed  value  of  the  universal  constant  q.  The 
temperature  variance  dissipation  rate  x  is  computed  directly  from  the  measurements 
and  e  is  treated  as  a  parameter  whose  value  is  specified  by  requiring  the  best  least 
squares  fit  of  (2.9)  to  the  measured  data.  This  method  based  on  the  Batchelor 
spectral  form  is  considered  accurate  within  a  factor  of  two.  I  have  attempted  to 
estimate  the  sensitivity /error  in  such  a  procedure  associated  with  using  either  the 
Batchelor  or  the  Kraichnan  model.  Several  subsets  of  our  data  set  (run  a240)  were 
selected  such  that  the  majority  of  the  points  either  belonged  to  pre-peak,  peak  or 
the  far  dissipation  range.  In  all  sets  the  dissipation  peak  was  resolved.  Using  earlier 
established  optimum  values  of  the  constant  q,  qKr-,  QBi  this  procedure  yielded 
the  relative  error  in  estimated  values  of  tjb  up  to  6%  when  using  the  Batchelor  form 
and  an  order  of  magnitude  less  in  the  Kraichnan  form  case.  Since  e  oc  7/^,  fitting 
the  Batchelor  form  to  experimental  data  can  account  for  at  least  25%  of  the  error  in 
c  estimates.  In  the  presence  of  systematic  measurement  noise  the  error  is  expected 
to  be  larger.  As  argued  by  Gibson  [23]  the  errors  in  estimating  the  dissipation  rates 
from  local  measurements  may  be  significantly  larger  if  intermittence  of  turbulence 
is  not  properly  taken  into  account.  Because  of  the  much  smaller  error  implied 
by  the  Kraichnan  formula  its  use  should  be  attempted  in  estimating  the  kinetic 
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energy  dissipation  rates  from  the  temperature  variance  measurements  in  oceanic 
environment. 

2.3.4  One  dimensional  scalar  spectra 

The  quantity  usually  measured  in  experiments  is  the  one  dimensional  scalar  spectrum 
Eie  which  is  related  to  its  three  dimensional  counterpart  E$  by  the  formula 

„  .  dEi$(k) 

^  =  -’‘Sk’  (2-27) 

and  normalized  as 

_  roo  l»oo 

Eu{k')dk'  =  jT  Ee{k)dk,  (2.28) 

We  can  obtain  the  one  dimensional  scalar  spectrum  derived  both  from  the  Batch¬ 
elor  and  from  the  Kraichnan  3-D  spectral  form  (2.9)  and  (2.11),  respectively 

[x{u^ll)hr]B\  ^  exp(-9(7/B/:)^)  -h  q^erf  {■qBk))  -  1).  (2.29) 

exp(-(6?)^BA:).  (2.30) 

In  Fig.  2.8  I  compare  our  calculated  one  dimensional  scalar  spectra  with  the 
above  formulas.  The  Batchelor  k~^  range  corresponds  to  the  flat  part  of  the  curves. 
Here  as  for  the  previous  three-dimensional  spectra,  Batchelor’s  prediction  agrees 
fairly  well  with  the  simulated  data  for  small  k  but  it  diverges  for  larger  wavenumbers 
{kqB  >1)5  while  Kraichnan ’s  model  describes  the  spectrum  very  well  in  the  entire 
range  of  simulated  wavenumbers. 

2.3.5  Location  and  universality  of  the  k~^  range 

The  remarkable  fact  is  that  Batchelor  scaling,  originally  associated  with  the  proper¬ 
ties  of  high  Reynolds  number  turbulence,  works  so  well  in  our  low  Reynolds  number 
simulations.  Batchelor  scaling  uses  the  rate  of  strain  of  the  Kolmogorov  eddies 
but  one  may  argue  that  in  our  forced  simulations  the  rate-of-strain  of  the 
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energy  containing  eddies  u' / Lp  is  more  appropriate.  However,  scaling  based  on  u' / Lp 
did  not  bring  the  data  to  a  self-similar  form.  The  success  of  the  Batchelor  scaling 
supports  the  arguments  put  forward  in  the  Introduction  where  I  propose  that  the 
Batchelor  range  spectrum  may  be  independent  of  Reynolds  number  because  of  the 
universality  of  the  dissipation  range  spectrum.  It  also  indirectly  supports  Gibson’s 
[22]  theory  which  predicts  that  the  Batchelor  scaling  should  apply  to  the  scalar  far 
dissipation  range  for  arbitrary  values  of  Pr,  not  necessarily  for  Pr  ^  1  as  first 
proposed  by  Batchelor  [3].  The  location  of  the  simulated  k~^  range  for  krjB  <  0.2, 
equivalent  to  roughly  krjK  <  1,  appears  to  be  at  variance  with  the  Batchelor  theory 
which  predicts  the  k  ^  range  for  wavenumbers  krjK  1.  However,  the  inconsistency 
is  removed  if  one  notes  that  the  velocity  scales  responsible  for  the  Batchelor  range 
behavior  should  be  those  that  provide  the  largest  rates-of-strain.  Since  the  spectral 
estimate  of  the  rate-of-strain  is  {E{k)ky/\  generally  the  dominant  scales  are  lo¬ 
cated  in  the  vicinity  of  the  maximum  dissipation  wavenumber  krjK  ss  0.2  rather  than 
around  the  Kolmogorov  wavenumber  krjK  «  1.  Moreover,  the  extent  and  the  loca¬ 
tion  of  the  Batchelor  range  is  consistent  with  the  following  Obukhov  type  argument 
[15].  The  scalar  variance  flux  F(k)  at  wavenumber  k  is  estimated  as 

F{k)  =  kFo{k)lT(k),  (2.31) 

where  r(A:)  is  the  turnover  time  for  scales  k.  In  Batchelor  theory  the  scalar  flux  and 
the  turnover  time  are  independent  of  k,  F{k)  =  x  and  r  ~  respectively, 

resulting  in 

(2.32) 

In  general  the  scalar  flux  is  expected  to  be  approximately  constant  in  the  range  of 
wavenumbers  with  negligible  scalar  variance  dissipation.  As  one  progresses  towards 
larger  wavenumbers  the  scalar  variance  flux  is  steadily  decreased  by  the  dissipation, 
with  the  largest  effect  around  the  scalar  dissipation  peak.  According  to  Fig.  2.6 
the  dissipation  peak  is  located  somewhere  between  (0.2  —  0.3)fc7/B.  Therefore,  the 
previous  arguments  suggest  that  the  Batchelor  range  should  be  observed  for  k  < 
0.2kT]B,  which  is  indeed  the  case.  In  general  the  experimental  results  are  consistent 
with  the  onset  of  the  Batchelor  range  at  around  krjK  —  0.05  for  Prandtl  numbers 
ranging  from  Pr  ~  1  to  Pr  =  100  [30].  In  our  data  the  beginning  of  the  k~^  range 
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can  be  observed  in  Fig.  2.9  where  the  scalar  spectra  multiplied  by  k  are  plotted 
after  normalizing  wavenumbers  with  the  Kolmogorov  length  •qx  rather  than  with 
the  Batchelor  length  qB-  With  this  normalization  I  obtain  a  family  of  self-similar 
curves,  one  for  each  different  value  of  the  Prandtl  number.  The  extent  of  the  k~^ 
range  increases  with  the  increasing  value  of  Pr  but  all  curves  converge  to  a  common 
origin  at  kqK  0.05. 

2.4  Conclusions 

We  have  examined  energy  and  passive  scalar  spectra  obtained  in  high  resolution  di¬ 
rect  numerical  simulations  at  three  different  Reynolds  numbers  and  for  three  Prandtl 
numbers  Pr  =  3, 5,  and  7. 

At  the  highest  Reynolds  number  Rx  =  77  the  energy  spectrum  has  a  short  inertial 
range  and  the  Kolmogorov  constant  is  consistent  with  experimental  data.  The  lower 
Reynolds  number  cases  do  not  exhibit  the  inertial  range  behavior.  The  maximum  in 
the  dissipation  spectrum  is  found  between  0.1  and  0.2kqK.  In  all  cases  the  spectra 
in  the  far  dissipation  range  for  kqic  >  0.2,  i.e.  beyond  the  dissipation  peak,  have  an 
universal  exponential  form  independent  of  Reynolds  number  and  with  parameters 
in  good  agreement  with  other  numerical  investigations.  The  universality  implies 
that  those  features  of  turbulence  that  axe  determined  by  scales  from  the  vicinity 
of  the  dissipation  peak  and  beyond  can  be  investigated  through  direct  numerical 
simulations  irrespective  of  the  value  of  the  Reynolds  number.  In  particular,  since 
the  Batchelor  spectrum  for  passive  scalars  with  Pr  >  1  is  assumed  to  be  controlled 
by  the  velocity  scales  from  the  vicinity  of  kqn  1,  this  suggests  that  passive  scalars 
with  Pr  >  1  advected  by  such  turbulent  velocity  fields  should  exhibit  universal 
behavior  consistent  with  the  predictions  of  Batchelor  [3]  and  Gibson  [22]. 

This  conclusion  is  confirmed  by  the  simulated  passive  scalar  spectra.  In  all  cases 
I  observe  the  appearance  of  the  Batchelor  k~^  range  followed  by  the  exponential 
decay  in  wavenumber  k  and  the  spectra  are  self-similar  under  the  Batchelor  scaling. 
However,  we  find  the  Batchelor  range  behavior  for  wavenumbers  kq^  <  1,  at  variance 
with  the  Batchelor  theory  which  predicts  the  k~^  range  for  wavenumbers  kqn  ^ 
1.  Our  results  suggest  that  the  velocity  scales  responsible  for  the  Batchelor  range 
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Figure  2.5:  Normalized  scalar  spectra  for  all 

all  Prandtl  numbers  (x-run  aal62,  o-run  a240,  *-run  xl62). 


runs 


Figure  2.7:  Normalized  scalar  dissipation  spectra  Eg{kT]B)[krig)^f[xii^f^Y^^ilB)  for 
run  a240  and  the  best  least  squares  fit  provided  the  Kraichnan  and  Batchelor  spec¬ 
tral  forms. 
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Figure  2.8:  Normalized  one  dimensional  scalar  spectra  Ei{kr^K){kr^K)l{x{t'leY''^'nB) 
for  all  runs  and  all  Prandtl  numbers  with  wavenumbers  nondimensionalized  by  the 
Kolmogorov  scale.  For  comparison  the  Kraichnan  (2.30)  formula  and  the  frequently 
used  Batchelor  expression  (2.29)  are  plotted. 
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Figure  2  9:  Normalized  scalar  spectra  E,(kr,K){kr,K)l(xWlc)'l^nB)  for  all  runs  and 
all  Frandtl  numbers  with  wavenumbers  nondimensionalized  by  the  Kolmogorov 
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behavior  reside  in  the  vicinity  of  the  maximum  dissipation  wavenumber  krjK  ^  0.2 
rather  than  around  the  Kolmogorov  wavenumber  krjK  ^  !• 

Comparisons  with  the  Batchelor  and  Kraichnan  functional  forms  of  the  scalar 
spectra  show  that  data  from  our  simulations  follow  closely  the  Kraichnan  form  while 
displaying  systematic  deviations  from  the  Batchelor  form.  In  the  context  of  oceano¬ 
graphic  procedures  which  estimate  the  kinetic  energy  dissipation  rates  e  from  the 
temperature  measurements  the  implication  is  that  the  Kraichnan  formula  will  pro¬ 
vide  more  accurate  estimates  than  the  Batchelor  formula. 

Our  numerical  results  contribute  to  the  ongoing  discussion  concerning  the  form 
and  universality  of  passive  scalar  spectra  in  turbulent  flows  for  Pr  >  1.  It  may  be 
useful  to  summarize  briefly  several  main  differences  between  our  assumptions  and 
experimental  conditions  and  those  of  other  studies  since  they  may  responsible  for  the 
different  conclusion  reached.  The  values  of  the  Prandtl  number  in  the  simulations 
(3,  5,  and  7)  are  greater  than  unity  but  much  less  than  those  of  some  experiments 
[47].  Similarly,  the  values  of  Reynolds  numbers  in  the  simulations  are  much  smaller 
than  in  many  experiments  but  comparable  to  others  [20].  The  simulations  are  very 
well  resolved  in  the  range  of  small  scales,  much  better  than  typical  experiments  or 
field  measurements  [53],  allowing  the  clear  distinction  between  the  Gaussian  and 
exponential  fall-off  at  large  k.  Contrary  to  other  numerical  simulations  [57]  which 
restrict  the  velocity  field  responsible  for  advecting  the  scalar  to  only  a  few  active 
modes,  I  use  the  full  velocity  field  obtained  from  the  Navier-Stokes  equations.  This 
may  enable  us  to  observe  the  k~^  range  at  smaller  resolutions  than  suggested  by 
other  numerical  models. 

Finally,  I  believe  that  analyses  of  the  velocity  and  scalar  field  using  physical 
space  representations  are  helpful  in  determining  the  precise  physical  mechanisms 
responsible  for  the  Batchelor  range  behavior  and  I  propose  that  such  analysis  be 
performed  in  the  future. 

2.5  Summary 

Direct  numerical  simulations  of  passive  scalars,  with  Prandtl  numbers  Pr  =  3,5, 
and  7,  advected  by  turbulence  at  three  low  Reynolds  numbers  were  performed.  The 
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energy  spectra  are  self-similar  under  the  R  I 

ior  consistent  with  many  other  investigatio  T""  '>^*>av. 

Reynolds  number  and  the  universal  exn  r.r  ‘he  highest 

numbers  in  the  dissipation  range.  I„  auT”  ''‘  °™  nil  Reynolds 

a  single  self-similar  curve  under  the  Satchel 

lowed  by  an  exponential  fall-off.  J  attribut  th  r*  ^“8^  «- 

for  our  low  Reynolds  number  flows  to  the  ‘‘“hRRy  of  the  Batchelor  scaling 

tra.  The  Batchelor  range  is  observed  Jw^Lr'' 
experimental  observations  but  smaller  tha  S®^®ral  agreement  with 

The  discrepancy  is  caused  by  the  fact  that'th'^  ^^l^sical  estimates, 

generation  of  the  Batchelor  range  are  in  .b  ‘he 

mum  energy  dissipation,  which  is  one  order  7'““^  "’“''“““her  of  the  max- 

wavenumber  used  in  the  classical  theory.  t1  77““'^' '^“'”°gorov 
scalar  spectra  proposed  by  Batchelor  and  K  ’  h  forms  of  passive 

results  and  I  found  that  the  Kraichnan  modT  ”  """  “““'ntions 

Batchelor  formula  displays  large  systemaf  d  """  ‘h'  The 

ntay  account  for  some  portion  of  the  obse  “d  this  fact 

calculations  which  estimate  the  kinetic  enir  oceanographic 

measurements.  <ijssipation  rates  from  temperature 
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Chapter  3 


A  comparison  of  the 

t-ne  eltects  on  near  forward 

»« pa«ic:::^ 

tne  ocean 


light 

in 


3.1  Introduction 

IS  changed  is  referred  to  scattering  i  individual  photon 

be  changed  in  the  near-for„,,<,  ,egil  ehheTr” 

geneities  or  particulates.  The  role  of  oa  f  i  ^  ^  interaction  with  water  inhomo- 
of  oceanic  waters  is  actively  studied  (' rs‘“ P'^P^rties 

Here  I  present  evidence,  obtained  from  1  ’  *“°™t. 

oratory  experiment,  that  turbulence  domirT'''  ^  Preliminary  lab- 

™‘ere.  --“ering  at  small  angles  in  oceanic 

Since  air  turbulence  limits  t  • 

optical  instruments  the  study  of  lighniT^  telescopes  and  other 

has  played  an  important  role  in  astronomi^^^  scattering  due  to  turbulence 

recent  progress  of  adaptive  optics,  which  can 

has  improved  the  resolution  of  earth  tel  atmospheric  turbulence, 

space  (turbulence  free)  counterparts.  A  re^^^^  ^  c  resolution  as  their 

atmosphere  is  presented  by  Kravtsov  [41]  ^  Propagation  in  the  turbulent 
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ect  of  turbuleoco  is  observed  al  smal)  „  ’  r  ’  ’  '  ‘hat  the 

«er,  the  .mportanee  of  turbulence  on  opfea  "™'’''  How- 

Th,s  controversy  is  in  large  part  due  to  th  controversial 

ogy  and  its  underwater  use  A  n!  of  the  instrument 

cause  computers  have  not  been  able  t  impeded  b 

addressed  than  for  the  atmosphere  I,  i,  acales  must  be 

i^‘-ons,  though  they  are  very  computer  out  numerical  simu 

“e  on  a  supercomputer  (5],  The  results  ZT"'  ““Po‘ation 

‘he  prevmus  chapter.  ‘hese  simulations  were  presented  in 

By  contrast,  the  role  of  particular  ■ 

received  much  attention.  This  is  becarse”  "  “sht  in  the  ocean  has 

on  are  responsible  for  the  variability  of  tir'^'r? P''J'‘“Plank- 

Z  *'>=-'>  «*ht,  thus  remoZT  “ot 

n  particulates  scatter  and  redistribute  1 J  r  “"®'  ""‘or  column, 

predominates  over  bacicscattering.  I„  f,,.  ®  ’““"riog  at  small  angles 

mattered  at  angles  less  than  3«  [631  Thl  -"“ored  light  “ 

rparched  subject  for  the  case  of  spherical  h^'*  O"  mattering  is  a  well 

scattering  on  particulates  is  guided  by  Mi  7“T  »‘udy 

The  importance  of  near-forward  scatter  "  '  ’ 

'  ®  on  underwater  imaging  th  W  by  Gordon  [9]  js  in 

;-0  p  surface,  and  on  ZZZZZZ  " 

Sma  l  angle  scattering  can  have  a  large  effect  'T'  °Ttical  sensors, 

ha^  become  increasingly  important  as  newTiLe  “t  ™s 

pable  of  imaging  to  live  optical  depths  and  ““‘rPments  are  ca- 

meters.  As  the  range  is  increased,  the  ‘^Tically  a  few  tens  of 

sensing  applications,  even  near  an  unscatt  d  7  '  '“‘o.  For  remote 

"  '  depend  on  small-angle  scatteri;  t^  77^"’  distributiol 

‘-■how-se,  optical  instruments  that  meas!  radiance 

pear  forward  scattering.  Additionally,  underwit™'"'"  ““  “f 

underwater  communication  is  degraded  by 
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pt.cal  pulse  stretching  which  depends  directly  on  It,,  1  f  , 

■ng.  The  transmission  rate  depends  on  th  ,  '  """ 

angle  scattering  processes.  '  amall 

I  propose  that  turbulence-induced  scaft.,-  ■  ■  c 

forward  scattering.  The  aim  of  this  cha  t  ^  for  most  near 

forward  scattering  on  partied  ^ 

^‘art  «rst  with  scattering  on  parluTalr^^""  »  -water.  I 

cattermg  on  particulates;  open  ocean  and 
coastal  example 

In  the  course  of  this  chapter  I  will  use  the  vol 

as  a  variable  to  describe  the  scatter]  a  f“notion  (VSF)  /J(a) 

describes  the  angular  distribution  of  th  ^  ^  medium.  The  VSF 

mtensity  df(a)  in  a  direction  o  per  nniUc  t 

irradianceBo  divided  by  the  incident 


I3(<x)  =  "W 
EodV' 


The  associated  scattering  coefficient  6  d  ’k  u 

and  scattered  from  the  incident  plan  '  removed 

ent  plane  wave  beam;  it  is  defined  as: 

/?(a)5m(a)dQr  ^2  2) 

It  should  be  noted  that  the  lieht  k« 

sorption  processes  and  the  corres  T  scattering  and  ab- 

P  -  <■  +  6,  (where  o  is  the  absorption  coefficient)  The  it  ?  ’  ^ 

rameteriaes  the  exponential  decay  of  the  fl  I  '  °  Pa- 

«M  -  '» 

■a  called  the  optical  thickness.  Even  thouvh  11  ^  "“■'^‘">™a,onal  quantity  cr 

wavelength-dependent,  for  the  purno  quantities,  a  and  6,  are  in  general 

length.  'will  consider  only  one  wave- 

The  VSF  and  6  are  considered  inherent  cmf  i 
fer  (56J,  these  are  properties  who,  ‘  Following  Preisendor- 

Perties  whose  magnitudes  depend  only  on  the  substances 
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g  t  Helds  that  may  pervade  it.  For  turbulence  this  becomes  less  clear  because  of  its 

1  column  can  only  be  characterized  by  a  given  turbulence 

for  s  ort  t, me  periods.  The  VSF  as  defined  above  provides  a  ood  description 
of  the  medmm  w.th  randomly  organized  structures  like  cells  where  we  can  assume 
hat  ax,al  symmetry  exists.  However,  it  fails  in  the  case  of  more  organized  media 
e  crystals  and  ,t  becomes  difficult  to  interpret  for  turbulent  media.  I  will  use 
nevertheless  because  .t  provides  a  useful  framework  to  compare  with  available 
theor.es  and  expenmental  evidence  of  scattering  in  the  ocean. 

VS  -herent  coefficients  and  in  the  case  of 

r  It  r” 7  ™bstances  are  present  in  the  medium,  the  global  VSF 

r  s  ts  from  the  sum  of  the  contributions  from  all  constituents  (as  long  as  the  singL 

scatter, ng  assumpt.on  is  fulfilled).  This  allows  me  to  calculate  a  total  VSF  in  the  case 

when  he  parfculate  ensemble  is  made  up  of  particles  of  various  size  categories  In 

general  scattering  by  spherical  particulates  is  exactly  described  by  Mie  theory.  The 
parameters  in  this  theory  are  refrartiv^a  •  j 

scatferma  f  11  i  Scattering  angle.  Fortunately, 

scattering  at  small  angles  can  be  ereatlv 

j;xr  . .  ,  r  ^  ^  mphfied  since  it  is  nearly  equivalent  to 

diifraction  by  apertures  of  equal  diameter  f^Ql  anri  :f  •  •  -r.* 

..  f  ^  meter  [59J  and  it  is  insensitive  to  variation  in 

e  refractive  mdex.  This  approximation  is  valid  within  a  couple  of  percent  no  to 

ang  es  o  0(10”).  The  typical  VSF  for  smaU  angles  due  to  diffraction  for  one  particle 
IS  shown  in  Figure  3.1.  p^iiuue 

For  our  compMsou  I  will  use  a  particle  size  distribution  that  is  characteristic 
s  a  wa  ers  [  6]  and  the  VSF  standard  for  open  ocean  clean  water  [32J.  The 
coital  particulate  distribution  was  measured  with  a  Coulter  counter  in  the  Bahama 

IT  -  7;  rT  of  -3  for  the  entire  (cu¬ 

mulative)  distribution.  I  then  calculated  the  particle  concentrations  in  individual 

Ttilr  t  ‘-0  corresponding  VSF  for  one 

!'  3.2.  The  resultant  VSF  for 

ang  es  up  to  0  la  shown  in  Figure  3.3.  This  VSF  exhibits  the  typical  properties  as 

experimen  al  y  obtained  by  Spinrad  (63),  namely  the  flat  plateau  for  angles  of  a  few 
egrees  and  decreasing  afterwards.  The  absolute  value  of  the  VSF  for  the  coastal 
ample  is  higher  by  about  a  factor  of  10,  as  would  be  expected  for  the  larger  concen¬ 
trations  and  the  presence  of  large  cells  and  is  consistent  with  the  lower  particulate 
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Figure  3.1:  Normalized  VSF  dnp  A-,ff 

sphere  of  size  parameters  i  =  20  and  i  ='lOo"  T  for  a  single 

the  near-forward  VSF  peak  is  given  by  A  =  fOO/xTn  . 


concentrations  used  bv  Soinrad  TKo  4.  •. 

VSF  at  small  angles  (a  0)  are  from  th  "  I 

^  ^  tile  largest  particulates  of  radius  1000//m. 

VSF  form  for  clear  open  ocean  waters 
proposed  by  Marten-Wells  [32].  This  VSF  has  h^  «  j  r 
the  choice  of  o,=  0.03re/  fo-mall  angles  and 


/?(a)  = 


27r(ag-f-a2)3/2- 


(3.3) 


3.3  Interaction  of  light  with  oceanic  temperature 
inhomogeneities 

Another  mechanism  of  light  scatterino-  ir,  fu  •  , 

•  I  .  .  ^  *11  t^e  ocean  is  the  scattering  on  turbulent 

:tr;rr “  ‘r.:r  ^  --  - 

f  f  ,  ®  ^  appear  (for  two  extreme  beam  diameters)  in 

two  forms:  for  wide  light  beams  (wider  than  small  temperature  scales)  the  ini  1  “ 

un,form  hght  beam  becomes  .granular-,  f.e,  exhibits  spatial  non-uniforl ^  d -  „ 

fo^mg  on  turbulent  inhomogeneities.  In  the  case  of  thin  light  beams  turb^Ll 

causes  beam  wandermg  or  .jitter’.  Due  to  technical  problems  it  is  experimentally 

cilfncult  to  measure  these  effects  Am/  i-n*  *  r  .  .  ^ 

likely  to  affect  thp  imperfection  in  the  optical  system  is 

we  hav  h  "measurements,  making  such  observations  unreliable  [74].  Thus 

we  have  chosen  to  quantify  the  efFerf  of  f  u  i  l  j  us 

simulation  of  th  f  k  i  n  urbulence  on  scattering  using  a  numerical 

mulation  of  the  turbulent  flow  field  and  of  light  propagation. 

3.4  Light  scattering  on  turbulent  flow 

sit  ("Jr  "T  "  con- 

c  ive  m  ex)  ,s  assumed  to  be  space  variant  and  the  effect  of  depolarization 
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3.2:  Size  distribution  of 
rdon’s  measurements  [19 


is  negligible.  If  backscatter  is  not  imoortanf  ft.  r  u. 

tKe  pa.aM.ed  (diff.iondi.e)  ^ 

o  / 

+  2f’nV  =  0,  (3 

Tf  l^ht’^ptopagtlf  v=l?; 

A  is  the  I  ft,  ’  ^  transverse  Laplacian  operator,  and  k  =  2n/X  where 

the  wave  ength  apd  .  the  deviation  of  the  local  .efcactive  index  f  Jts  ll 

however  I  haveTo^Z' equlTio^^ 

Further  •  1-fi  •  ^  ^  developed  numerical  algorithm, 

propagation  over  short  distanc  s  hi  I  t  “ 

,  lations  0.25  „.)  obeys  («)-/C  „  1  "  <“  <>-  -Icu- 

— eofthe«o.,ypicallyof0l’^^^^^^^^^^^^^^ 

‘-“-r  “  ,.L  :ii ; 

V’(X, y)  =  .4„(s:, y) exp  (^-ik £  y, 

where  A,  ,s  mitial  value  of  i,(x,  y,  a)  at  a  =  0.  This  solution  constitutes  the  first 

negl  ”  Z  Si-  ‘ire  second  derivatives  have  been 

n^feted  n  th,  approx.nrat.on  only  changes  in  phase  can  be  calculated.  This  is 
a  sufficient  simplification  for  the  oresent  •  t 

light  scattering  r.  1  f  ^  Since  I  am  mostly  interested  in  the 

g  scatter  ng  angle  for  propagation  distances  over  a  few  tenths  of  cm.  To  verify  if 
he  geometncaJ  optics  approximation  (E,.  3.5,  used  in  this  work  is  appropril  H 
mvest, gated  parameter  range  of  distances  and  turbulent  flows  we  have  validated 

ThlTirl  Eq.  (3.4)  [61. 

a  ig  p  asc  at  the  exit  of  the  computational  volume  can  be  obtained  with 

ins°r  “ti  with  knowledge  of  the  spatial  distribution  of  the  refractive  index 

ns  de  he  investigated  volume.  Since  there  is  a  simple  relation  between  tempZt! 

and  refractive  index  the  spatial  distribution  of  the  temperature  field  is  requZ 
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I!"™”'"' «o„ 

number  of  the 

are  relevant  for  light  scattering  For  exa  T”“'  ''  ^ow 

‘ure  variability  on  scales  ranging  fro„  0(l7l^)  7  ““““ 

approach  to  determine  the  optically-re  Jant  r  r‘“  ' 

geometrical  optics.  It  can  be  shown  that  th  7"*'  “  obtained  from 

-i-e  (ot),  betwXtutT 

ray  direction)  per  unit  distance  dr  depend  ,  “■‘-I 

temperature  distribution  in  the  absence  of  th”  I  ““ 

the  flow  i.e.;  o  arge  scale  temperature  gradients  in 

d(a^)  roo 

~  dr  ~  ^  L  ^Tkdk 

where  dr  is  a  smail  increment  along  the  linbt  b  • 

propagation  coordinate,  and  EM)  is  the  t  w  P^Pon^ioular  to  the  light 

■k)  of  the  temperature  field  Ttr)  (see  tb  '  (in  Fourier  space 

tion  to  the  integral  on  the  right  hand  sid  chapter).  The  largest  contribu- 

the  observed  oceanic  range  of  c  val  propagation  distances  of  0(m)  and 

the  turbulent  field  (5J.  Therefore  any  7um7r”7l  ^""‘7  “ 

tive  index,  here  due  to  the  temperature  fl  ''frac- 

correctly  to  account  for  the  most  important  77““’ 

We  chose  a  model  of  ho„  '”P‘’^‘“t  scattering  events. 

“homogeneities  of  the  temperlturrfieirXl”  Sow  to  describe  the 

which  there  exist  theoretical  solutions  H  '  ®™Piost  turbulence  model  for 

stable  temperature  gradients  which  mav  b””'""’  “  '"S' 

boundary  layer  flows.  ^  Present  m  actual  ocean  flows  such  as 

“  77,^17:7“’  '‘“““‘-■'■s 

at  small  angles  on  RII 

As  demonstrated  earlier  from  a  Ueht 

to  scattering  comes  from  the  smallesn7l77f“tr'r‘'‘'”‘“‘’  “““^utiou 

e  emperature  Bt  distribution  (Eq. 
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temperature  ae  a  passive  '‘‘’7  '  "ir  »"  ““taining  the 

the  temperature  variance  17' t"'  T'*  ^  distribution  of 

the  Kraichnan  model,  E,.  (^SoJ 

whic7r 01?;^;^ :  ^T  ^  ^ 

kinetic  energy,  e,  determines  the  spatial  size  ,, 

temperature  structure  THp  nK  •  i  ^7  (  /e)  ,  of  the  smallest 

ture.  ihe  physical  meaning  of  y  (y  =  p  u2ju\  ■ 

at  which  the  temperature  flnrf.,  r  ^  ^  ° 

perature  fluctuations  present  in  the  flow  dissioate  f-?!!  v 

optical  perspective,  ^  has  the  following  interpretation-  la  t  “ 

to  stronger  temperature  vraH'  r  'n'nrpretation.  larger  values  correspond 

flow.  This  situation  is  oi  smallest  temperature  structure  of  the 

optical  properties  of  a  turLl7t  h  summarize,  the 

by  two  parameters:  ;r  which  ezpre  T”'”"'  characterized 

e  which  is  inversely  propo  T  T 

ranges  from  10-  «  CVs  Tf  w  7  77'’' 

the  mid  water  column  (16)  For  e  th  T  “ 

lioj.  for  e  the  range  is  from  10~‘*m2/<!3  _  r,-  i 

upper  layer  to  10““m2/e3  ^  fairly  energetic 

y  r  lo  lu  m  /s  in  mid  water  column  [1] 

or  rircoTuii  is  :“?u::rh^  ^  — 

plot  of  observed  y  us.  e  is  sh  •  Z 

us.  r  are  fairly  uncorrelated.  k 

3.7  Numerical  experiment  of  light  scattering  on 
turbulent  flow 

As  established  in  the  previous  chanter  th  t 

scales  are  independent  of  the  th  '  '  ‘™Pcrature  variance  spectra  at  small 

follows.  First  I  scale  r?  n  -  - 

2.30  given  the  values  of  y^Tc  aid  ^• 

the  refractive  index  n'(x  u  z)  Th  I  '1^  distribution  of 

,y,  ).  Then  I  illuminate  the  computational  volume  with 
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turbulent  blob -£' 


Figure  3.4: 


Optical  interpretation  of  the  tnrbnlent 


parameters  and  e. 
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Figure  3.5:  Profiles  of  e,  x,  temperature  and  salinity  from  the  Oregon  coast  (courtesy 
of  J.  Mourn). 


a  plane  parallel  wave  by  assigning  a  constant  value  to  the  function  tf}(x,y,z  =  0). 
Since  I  illuminate  our  computational  volume  (240^  size)  with  a  plane  parallel  wave, 
the  phase  and  amplitude  of  the  light  is  the  same  throughout  the  entrance  plane.  We 
then  use  Eq.  (3.5)  (or  Eq.  (3.4)  for  verification)  to  estimate  the  resultant  phase 
distortion  0(x,t/)  of  light  on  the  exit  side  of  the  computational  volume. 

The  imaginary  value  taken  by  ^  at  each  x,  y  position  describes  how  the  initially 
plane-parallel  wavefront  has  been  distorted  by  passage  through  the  temperature 
inhomogeneities  of  the  computational  volume.  The  real  part  (different  from  Ao  only 
if  I  solve  Eq.  (3.4))  describes  the  amplitude  fluctuations.  I  show  an  example  of  the 
calculated  distribution  of  light  intensity  at  the  exit  estimated  using  Eq.  (3.4)  in  [6]. 
It  can  be  seen  there  that  the  initially  homogeneous  light  beam  becomes  spatially 
inhomogeneous  due  to  the  effect  of  turbulence-induced  scattering  alone. 

The  imaginary  part  of  ij}{x,  y)  carries  information  about  scattering  angle  of  the 
light  beam.  The  vector  normal  to  the  isosurface  z  =  rp{^x,y)  is  N'. 


N  = 


dij)  dxf) 
dx'  dy' 


(3.7) 


and  is  parallel  to  the  outgoing  ray  direction.  Since  the  exit  surface  has  240^  different 
elements  (of  different  N),  our  resolution  is  limited  to  this  number  when  determining 
the  VSF.  The  final  calculation  of  the  VSF  is  done  in  a  straightforward  manner  from 
Eq.  3.1.  The  calculated  VSF  for  extremal  (large  x)  conditions  of  x  and  e  in  the 
coastal  ocean  is  given  in  Figure  3.7.  From  Figure  3.7  it  can  be  seen  that  the  VSF 

for  any  given  x  and  e  can  be  approximated  by  a  Gaussian  function,  determined  by 
flow-dependent  parameters: 


X,  a)  =  Vo(e, x)  ■  exp(-(a/ao(€, x))^),  (3.8) 

where  K,  describes  the  maximum  value  of  VSF  at  a  =  0  and  ao  is  the  half-width. 
The  scattering  angle  variance  ((a^))^/^,  calculated  with  respect  to  its  initial  di¬ 
rection,  is  dlso  estimated  as  a  diagnostic  quantity.  The  value  of  the  scattering  angle 
variance  depends  on  the  propagation  distance  so  I  use  a  constant  optical  path  0.25m. 


47 


Figure  37:  Simulated  extreme  VSF  for  a  range  of  turbulent  flo« 

egime  or  sea  surface:  a  range  of  y  and  constant  t  =  10-‘»  mVs^ 
denotes  particulate  VSF.  rn  js  . 


in  the  coastal 
The  solid  lines 


■18 


To  obtain  scattering  angle  variance  ((a=))‘/2  fo,  longer  propagation  distances  i,  the 
following  formula  can  be  used  [45]: 


X.  £)  =  CoMx,  e)  ■ 


(3.9) 


where  Lq  in  case  of  my  simulations 
scattering  angle  variance  has  been 
shown  in  Table  3.1. 


IS  0.25  m.  The  values  of  Vo,  qq  and  Co.25(x,  e) 
calculated  for  different  flow  conditions  and  is 


49 


Table  3.1: 
length. 


wmuKM 

m 

ao  [deg] 

^0.25 (y,^)  fradl 

-10 

-2 

2*103 

1.1  *  10-* 

10-3 

-10 

-4 

2*105 

1.1  *  10-2 

10-“ 

-10 

-6 

2*10^ 

1.1  *  10-3 

10-5 

-8 

-2 

2  *  10® 

5  *  10-2 

1 

o 

* 

-8 

-4 

2*10’’ 

5  *  10-3 

5  *  10-5 

-8 

-6 

2*  10^ 

5  *  10-^ 

5  *  10-5 

-6 

-2 

2*  10® 

1.1  *10-2 

2.5  *  10-“ 

-6 

-4 

2*10* 

1.1  *  10-3 

2.5  *  10-5 

-6 

-6 

2*  10*® 

1.1  *  lO-'* 

2.5  *  10-5 

VSF(K„,ao)  and  the  scattering  angle  variance  ({a^))‘/2  for  0.25  m  path- 
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3.8  Comparison  of  VSF  for  particulates 
and  turbulence  in  the  coastal  zone 

The  VSF  corresponding  to  particulates  obtained  earlier  can  be  compared  with  the 
VSF  corresponding  to  energetic  turbulence  of  the  coastal  ocean.  From  Table  3.1  we 
can  see  that  the  turbulent  flow  scatters  light  most  when  x  is  large  and  c  is  low.  If 
we  assume  that  x  is  driven  by  the  mean  gradient  and  a  turbulent  time  scale  (u'/T), 
where  u'  is  the  typical  rms  current  velocity  and  T  is  a  characteristic  length  scale 
(for  example  the  depth  of  the  water  column  or  a  mixed  layer  depth  is  of  0(10  m)  in 
the  coastal  region),  we  obtain  the  relationship  x  ST^iu'/L),  where  8T  is  a  mean 
temperature  gradient.  From  this  argument,  I  obtain  values  for  x  of  IQ-^  and 

lower.  This  is  consistent  with  measurements  made  by  Farmer  and  Gemmrich  [19]  in 
the  surface  layer  (upper  few  meters)  of  open  ocean.  The  coastal  value  of  c  can  be 

given  analogously  by  e  ~  u^/L  [70]  and  can  be  as  high  as  10-^mVs^  and  as  low  as 
in  the  presented  data  set  shown  in  Fig.  3.5, 

Using  the  higher  limit  values  of  x,  I  have  estimated  the  VSF  for  a  range  of  x  and 
for  constant  e  —  10  m  /s  ,  value  which  is  likely  to  be  encountered  in  the  coastal 
region  or  in  the  open  ocean  (Figure  3.7).  For  comparison  the  VSF  corresponding  to 
the  particulates  of  the  coastal  ocean  and  open  ocean  are  also  plotted.  It  is  clearly 
shown  that  in  our  simulations  of  the  coastal  zone,  turbulence  scattering  can  dominate 
to  angles  as  large  as  5  .  The  implications  for  underwater  visibility  applications  are 
considerable  as  the  scattering  angle  variance  is  up  to  0(0.001  rad)  (Table 

3.1)  implying  that  over  a  range  of  10m  the  smallest  detail  which  can  be  resolved 
has  a  size  of  O(O.lm).  Since  the  surface  layer  of  the  open  ocean  shares  turbulence 
characteristics  with  those  of  the  coastal  zone,  the  VSF  also  describes  the  surface  open 
ocean.  Turbulence  dominates  scattering  for  angles  of  10^  The  visibility  through  the 
turbulent  surface  of  thickness  10m  is  also  going  to  be  limited  (in  absence  of  other 
effects)  to  detail  the  extent  of  0(0.  Im). 
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3.9  Experimental  validation 


To  qualitatively  validate  the  calculated  VSF.  I  carried  out  a  simple  laboratory  ex- 
penmen  at  the  laboratory  facility  of  WetLabs  (Philomath,  Oregon).  The  source 

th  b  h  T'  ir  optical  measurements  were  made  using 

wltllf  (VABAM)  developed  by 

Turbulence  for  the  experiment  was  generated  by  a  simple  convective  cell  created 
a  ea  e  e  ement  on  the  bottom  and  a  cooling  element  on  the  top.  The  applied 

hat  the  flow  was  fully  turbulent.  The  corresponding  ^  ranged  from  10-  down  to 
/s  and  r  values  were  between  10-  to  lO-mVs’.  These  wide  ranges  result 
rom  e  temperature  measurement  system  and  the  poor  temporal  and  spatial  uni- 
form.ty  of  the  heatmg/cooling  setup.  The  large  range  also  implies  that  the  turbulent 
as  m  omogeneous  and  anisotropic  and  therefore  not  well  controlled, 
e  light  path  m  VABAM  starts  at  the  flash  lamp  source,  goes  through  a 
monochromator  (.  420nm)  and  then  a  0.7  mm  pinhole  aperture  to  be  collimated 

vol  ^  Th^h  ■  c**’™ 

e  earn  is  ocuMed  by  another  200  mm  lens  (with  iris  aperture  is  driven 
y  a  s  epper  motor)  and  impinges  upon  a  diffuser  plate  followed  by  a  large  area 
detector.  The  minimum  aperture  diameter  of  the  iris  of  the  receiver  lens  defines  the 
acceptance  angle  of  the  receiver.  It  is  1.2  mm,  which  results  in  a  0.12«  acceptance 

angle  in  water.  The  step  resolution  of  the  stepper  motor  is  approximately  1  mm; 

thus  the  resolution  of  the  ^  i  r 

ptance  angle  ranges  from  approximately  0.1  to  3® 

1  ne  sampling  frequency  is  1  Hz. 

The  presence  of  turbulence  in  the  sample  volume  caused  very  intermittent  read- 
mgs  at  all  apertures;  the  largest  variability  was  observed  for  the  smallest  aperture 
The  intermittency  of  the  detected  light  is  shown  in  Figure  3.8.  Because  of  this  ex¬ 
treme  intermittency,  a  statistically  stable  value  of  the  VSF  could  only  be  obtained 
by  averaging;  We  calculated  the  VSF  using  data  from  a  run  of  over  an  hour  time 

perio  .  e  average  becomes  stable  after  averaging  hundreds  of  readings  over  a 
couple  of  minutes  and  is  shown  in  Figure  3.9. 


52 


Relative 
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Figure  3.8;  Time  series  of  VABAM  detector  readings  for  two  apertures.  The  smaller 
aperture  corresponds  to  the  larger  variance 
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Figure  3.9:  Experimentally  measured  VSF  (circles)  in  the  presence  of  turbulence 
temperature  inhomogeneities  and  the  simulated  VSF  (stars)  for  the  corresponding 
values,  of  e  and  x- 
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To  evaluate  the  effect  of  particles,  15  fim  beads  were  added  to  the  fairly  weak 
turbulent  flow.  In  another  experiment  this  concentration  resulted  in  a  visible  increase 
of  the  VSF  (consistent  with  Mie  theory)  in  runs  with  no  turbulence  [50].  However, 
in  the  weak  turbulent  regime  the  VSF  remained  larger  than  Mie  predictions  and 
there  was  no  change  in  VSF  after  adding  the  beads. 

It  should  be  kept  in  mind  that  in  theory  the  VSF  is  obtained  as  a  limit  when 
the  pathlength  goes  to  zero.  Usually  measurements  are  done  over  a  small  volume  V 
and  the  outcome  of  the  calculated  VSF  should  not  depend  on  the  volume  size.  In 
the  case  of  the  scattering  on  turbulence  this  is  not  true.  Since  turbulence  has  an  as¬ 
sociated  length-scale  77,  every  near-forward  VSF  measurement  will  depend  not  only 
on  turbulent  quantities  tj  and  x  but  also  on  the  size  of  the  volume  itself.  This  situ¬ 
ation  makes  comparison  of  near- forward  VSF  measurements  from  different  volumes 
difficult,  especially  when  made  from  different  instruments. 

Even  though  our  experimentally  measured  VSF  agrees  qualitatively  with  the  sim¬ 
ulation  it  is  difficult  to  compare  the  two  for  several  reasons.  Turbulence  simulations 
correspond  to  homogeneous  isotropic  turbulence,  the  values  of  x  ^  2ire  precisely 
determined,  and  I  have  only  one  realization  of  the  temperature  field.  By  contreist, 
the  turbulence  of  the  experiment  is  not  homogeneous;  I  cannot  provide  precise  values 
of  X  or  €,  and  many  realizations  of  the  temperature  field  are  averaged  throughout 
the  course  of  the  experiment  into  one  VSF  value.  We  know  from  the  simulations 
that  the  slope  of  VSF  is  very  sensitive  to  the  value  of  the  turbulent  parameters  x 
and  e  which  precludes  exact  predictions  for  our  experiment. 

The  most  important  outcome  of  the  experiment  is  the  confirmation  that  turbu¬ 
lence  indeed  dominates  scattering  at  small  angles  of  0(1®)  and  that  scattering  can 
be  well  modeled  using  our  numerical  DNS  databaise.  A  striking  result  is  that  the 
value  of  the  VSF  at  small  angles  is  highly  intermittent.  Any  meaisurement  technique 
which  relies  on  angles  of  0(1®)  or  less  must  develop  an  appropriate  averaging  scheme. 
This  intermittency  can  only  be  attributed  to  turbulence.  Another  major  conclusion 
was  that  turbulent  scattering  was  found  to  be  strong  enough  to  account  for  the  ob¬ 
served  values  of  the  near-forward  VSF.  The  example  of  in  situ  measurement  affected 
by  near-forward  light  scattering  on  turbulence  is  the  measurement  of  beamC.  To 
minimize  influence  of  turbulence  the  several  instantaneous  beamC  readings  should 
be  averaged  into  one  value. 
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In  spite  of  the  large  uncertainties  of  the  experimental  turbulence  quantities,  these 
values  are  highly  consistent  with  the  numerically  simulated  VSF.  This  experimental 
confirmation  of  the  predicted  VSF  can  help  us  to  understand  light  propagation  for 
example  in  the  coastal  zone. 

3.10  Conclusions 

•  Turbulent  scattering  dominates  particulate  scattering  up  to  5°  in  the  coastal 
zone  and  to  10®  in  the  upper  layer  of  the  open  ocean. 

•  The  turbulence  can  limit  underwater  visibility  to  details  of  size  0(0. Im)  for  a 
range  of  10  meters. 

•  Since  turbulent  flows  axe  highly  intermittent  this  results  m  high  variance  in 
the  value  of  the  VSF  at  small  angles. 

•  Measurement  of  the  VSF  in  the  near-forward  direction  depends  on  the  volume 
and  flow  characteristics  of  the  sampling  volume. 

•  Our  numerical  simulations  were  supported  by  the  preliminary  experiment.  A 
future  experiment  should  include  a  more  comprehensive  control  of  turbulent 

parameters. 

3.11  Summary 

We  have  investigated  numerically  and  experimentally  the  scattering  of  a  light  beam 
at  small  angles.  I  compare  the  effects  of  scattering  by  typical  coastal  and  oceanic 
particulate  distributions  with  turbulence-induced  scattering.  I  find  that  turbulence 
dominates  scattering  up  to  5-10®  for  the  coastal  region  and  for  the  upper  layer  of  the 
open  ocean.  Our  preliminary  experiment  confirms  predictions  that  the  frequently 
observed  high  in  situ  values  of  the  volume  scattering  function  at  small  angles  are 
related  to  scattering  on  turbulent  inhomogeneities  in  seawater.  Furthermore,  the 
initially  homogeneous  light  becomes  spatially  inhomogeneous  and  highly  intermit¬ 
tent  after  propagating  through  the  turbulent  region.  This  intermittence  will  affect 
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single  realization  measurements  by  increasing  the  variance.  In  addition  I  confirmed 
that  the  mean  scattering  angle  depends  strongly  on  x- 
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